rm(list = ls())
library(pacman)
pacman::p_load(readstata13,foreign,dplyr,data.table,lfe,ggplot2,wesanderson,biostat3)

d <- read.dta13('data_final.dta')

# Regression function
####################################################################################################

gen.reg <- function(outcome, ctrl, spec, df){
  if(outcome=='Knowledge'){fmla <- paste(outcome, '~ T')}
  if(outcome=='Behavior'){fmla <- paste(outcome, '~ T*List')}
  if(ctrl==1){fmla <- paste(fmla, '+ qualtrics + urban + female')}
  if(spec==1){fmla <- paste(fmla, '| block + Week | 0 | Week_id')}
  if(spec==2){fmla <- paste(fmla, '| id + Week | 0 | Week_id')}
  
  fmla <- as.formula(fmla)
  out <- felm(fmla, data=df)
  
  if(outcome=='Knowledge'){
    coef <- out$coefficients[rownames(out$coefficients)%in%c('T','T:List')]
    se <- out$cse[rownames(out$coefficients)%in%c('T','T:List')]
  }
  
  if(outcome=='Behavior'){
    coef <- out$coefficients[rownames(out$coefficients)%in%c('List','T:List')]
    se <- out$cse[rownames(out$coefficients)%in%c('List','T:List')]
    coef <- c(coef[1], coef[1]+coef[2])
    lin <- lincom(out,c("List+T:List"),eform=F)
    se <- unlist(c(se[1], as.vector((lin[3]-coef[2])/1.96))) # get SEs for the linear combination
  }
  
  return(cbind(coef=coef,se=se,spec,ctrl))
}


# Knowledge
####################################################################################################

out_1_0 <- as.data.frame(gen.reg('Knowledge',0,1,d))
out_1_1 <- as.data.frame(gen.reg('Knowledge',1,1,d))
out_2_0 <- as.data.frame(gen.reg('Knowledge',0,2,d))
out_2_1 <- as.data.frame(gen.reg('Knowledge',1,2,d))

out_k <- dplyr::bind_rows(out_1_0, out_1_1, out_2_0, out_2_1)
out_k$ctrl <- ifelse(out_k$ctrl==0,'Without\ncontrols','With\ncontrols')
out_k$ctrl <- as.factor(out_k$ctrl)
out_k$ctrl <- relevel(out_k$ctrl, ref='Without\ncontrols')
out_k$spec <- ifelse(out_k$spec==1,'Randomization Block Fixed Effects','Broadcast List Fixed Effects')
out_k$spec <- as.factor(out_k$spec)
out_k$spec <- relevel(out_k$spec, ref='Randomization Block Fixed Effects')

make.plot <- ggplot(data=out_k, aes(x=ctrl))+
  theme_bw()+xlab(NULL)+ylab('Treatment effect')+
  facet_grid(cols=vars(spec))+
  geom_col(aes(y=coef),colour='gray50',fill='gray80',width=0.75)+
  geom_point(aes(y=coef),pch=3,size=2)+
  geom_errorbar(aes(ymin=coef-(1.96*se),ymax=coef+(1.96*se)), width=0.2)+
  geom_hline(aes(yintercept=0))+
  theme(strip.background =element_rect(fill="white"))

pdf('Figures/knowledge.pdf',width=5,height=3)
print(make.plot); dev.off()
ggsave("Figures/knowledge.eps", width = 5, height = 3, units = "in")

out_1_female <- as.data.frame(gen.reg('Knowledge',0,1,d[d$female==1,]))
out_1_male <- as.data.frame(gen.reg('Knowledge',0,1,d[d$female==0,]))
out_1_urban <- as.data.frame(gen.reg('Knowledge',0,1,d[d$urban==1,]))
out_1_rural <- as.data.frame(gen.reg('Knowledge',0,1,d[d$urban==0,]))
out_1_week1 <- as.data.frame(gen.reg('Knowledge',0,1,d[d$Week==1,]))
out_1_week2 <- as.data.frame(gen.reg('Knowledge',0,1,d[d$Week==2,]))

out_k2 <- dplyr::bind_rows(out_1_female, out_1_male, out_1_urban, out_1_rural, out_1_week1, out_1_week2)
out_k2$spec <- c('Female','Male','Urban','Rural','Week 1','Week 2')

make.plot <- ggplot(data=out_k2, aes(x=spec))+
  theme_bw()+xlab(NULL)+ylab('Treatment effect')+
  geom_col(aes(y=coef),colour='gray50',fill='gray80',width=0.75)+
  geom_point(aes(y=coef),pch=3,size=2)+
  geom_errorbar(aes(ymin=coef-(1.96*se),ymax=coef+(1.96*se)), width=0.2)+
  geom_hline(aes(yintercept=0))+
  theme(strip.background =element_rect(fill="white"))

pdf('Figures/knowledge2.pdf',width=5,height=3)
print(make.plot); dev.off()
ggsave("Figures/knowledge2.eps", width = 5, height = 3, units = "in")


# Behavior
####################################################################################################

out_1_0 <- as.data.frame(gen.reg('Behavior',0,1,d))
out_1_1 <- as.data.frame(gen.reg('Behavior',1,1,d))
out_2_0 <- as.data.frame(gen.reg('Behavior',0,2,d))
out_2_1 <- as.data.frame(gen.reg('Behavior',1,2,d))

out_b <- dplyr::bind_rows(out_1_0, out_1_1, out_2_0, out_2_1)
out_b$ctrl <- ifelse(out_b$ctrl==0,'Without\ncontrols','With\ncontrols')
out_b$ctrl <- as.factor(out_b$ctrl)
out_b$ctrl <- relevel(out_b$ctrl, ref='Without\ncontrols')
out_b$spec <- ifelse(out_b$spec==1,'Randomization Block\nFixed Effects','Broadcast List\nFixed Effects')
out_b$spec <- as.factor(out_b$spec)
out_b$spec <- relevel(out_b$spec, ref='Randomization Block\nFixed Effects')
out_b$Condition <- rep(c('Control','Treatment'),4)

make.plot <- ggplot(data=out_b, aes(x=ctrl))+
  theme_bw()+xlab(NULL)+ylab('List experiment treatment effect')+
  facet_grid(cols=vars(spec))+
  geom_col(aes(y=coef,color=Condition),fill='gray80',width=0.6, position = position_dodge(width = 0.75))+
  geom_point(aes(y=coef),pch=3,size=2, position = position_dodge2(width = 0.75))+
  geom_errorbar(aes(ymin=coef-(1.96*se), ymax=coef+(1.96*se)),width=0.75, 
                position = position_dodge2(padding=0.5))+ #width=0.5,
  theme(strip.background =element_rect(fill="white"))+
  geom_hline(aes(yintercept=0))+
  scale_color_manual(values=wes_palette(n=2, name="Royal1"))

pdf('Figures/behavior.pdf',width=5,height=3)
print(make.plot); dev.off()
ggsave("Figures/behavior.eps", width = 5, height = 3, units = "in")

out_1_female <- as.data.frame(gen.reg('Behavior',0,1,d[d$female==1,]))
out_1_male <- as.data.frame(gen.reg('Behavior',0,1,d[d$female==0,]))
out_1_urban <- as.data.frame(gen.reg('Behavior',0,1,d[d$urban==1,]))
out_1_rural <- as.data.frame(gen.reg('Behavior',0,1,d[d$urban==0,]))
out_1_week1 <- as.data.frame(gen.reg('Behavior',0,1,d[d$Week==1,]))
out_1_week2 <- as.data.frame(gen.reg('Behavior',0,1,d[d$Week==2,]))

out_b2 <- dplyr::bind_rows(out_1_female, out_1_male, out_1_urban, out_1_rural, out_1_week1, out_1_week2)
out_b2$spec <- c('Female','Female','Male','Male',
                 'Urban','Urban','Rural','Rural',
                 'Week 1','Week 1','Week 2','Week 2')
out_b2$Condition <- rep(c('Control','Treatment'),6)

make.plot <- ggplot(data=out_b2, aes(x=spec))+
  theme_bw()+xlab(NULL)+ylab('List experiment treatment effect')+
  geom_col(aes(y=coef,color=Condition),fill='gray80',width=0.6, position = position_dodge(width = 0.75))+
  geom_point(aes(y=coef),pch=3,size=2, position = position_dodge2(width = 0.75))+
  geom_errorbar(aes(ymin=coef-(1.96*se), ymax=coef+(1.96*se)),width=0.75, 
                position = position_dodge2(padding=0.5))+ #width=0.5,
  theme(strip.background =element_rect(fill="white"))+
  geom_hline(aes(yintercept=0))+
  scale_color_manual(values=wes_palette(n=2, name="Royal1"))

pdf('Figures/behavior2.pdf',width=5,height=3)
print(make.plot); dev.off()
ggsave("Figures/behavior2.eps", width = 5, height = 3, units = "in")


# Trust in different media sources
####################################################################################################

trust_dat <- read.csv("trust_sources.csv", sep = ",")
vals <- colMeans(trust_dat[,2:9], na.rm = FALSE, dims = 1)

rev_vals <- rev(vals)

names <- c("From family and friends", 
           "From the Government",
           "From international organizations and agencies",
           "From NGOs, or CSOs",
           "Mention a news source",
           "Mention a doctor as a source",
           "Mention a government source",
           "None of these")

rev_names <- rev(names)

par(mai=c(1,4,1,1))
barplot(rev_vals,
        main = "Trust in each source (select up to three)",
        xlab = "Percentage",
        ylab = "",
        names.arg = rev_names,
        col = "grey",
        horiz = TRUE,
        xlim = c(0, 0.8),
        las = 1)
